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Abstract We describe the numerical modeling of current 
flow in graphene heteroj unctions, within the Keldysh Lan- 
dauer Non-equilibrium Green's function (NEGF) formalism. 
By implementing a ^- space approach along the transverse 
modes, coupled with partial matrix inversion using the Re- 
cursive Green's function Algorithm (RGFA), we can sim- 
ulate on an atomistic scale current flow across devices ap- 
proaching experimental dimensions. We use the numerical 
platform to deconstruct current flow in graphene, compare 
with experimental results on conductance, conductivity and 
quantum Hall, and deconstruct the physics of electron 'op- 
tics' and pseudospintronics in graphene p — n junctions. We 
also demonstrate how to impose exact open boundary condi- 
tions along the edges to minimize spurious edge reflections. 

Keywords NEGF, RGFA, graphene, electron 'optics' 
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1 Introduction 

Ever since its inception 1 1 1, graphene as a material has occu- 
pied a unique place between one dimensional pi-conjugated 
organic molecules and three dimensional bulk crystalline 
solids with well defined bandstructures. The quasi-ballistic 
scattering lengths ||2j|3|, photon like dispersion |4 | and chi- 
ral electron flow | 5 1 together promote nontrivial transport 
physics in graphene, opening up various device possibilities 
that necessitate detailed numerical modeling. While the pho- 
tonic bandstructures argue for a continuum model of mod- 
est sized device segments, scattering at heteroj unction inter- 
faces |[6 1 and edges requires a detailed atomistic treatment. 
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The aim of this paper is to describe practical techniques for 
simulating quantum transport through graphene based het- 
eroj unctions, compare various transport regimes with exper- 
iments associated with the chiral nature of electron transmis- 
sion at graphene heteroj unctions. We also point out some of 
the unique electronic properties of bilayer graphene, but un- 
less otherwise stated, the results are for single layer graphene. 



2 Simulation platform 

The Non-equilibrium Green's function (NEGF) formalism 
|7| provides a unified, 'bottom-up' platform for modeling 
quantum flow of electrons. Indeed, it has been successfully 
applied to understand transport physics in materials and sys- 
tems as diverse as organic molecules |8|, carbon nanotubes 
|[91[TQ|, graphene |^-^16J, silicon nanowires |T7|j22| , spin- 
tronics, nanomagnets ||23j|24j, nanoscale phonon transport 
|25 -27|. A challenge however is the sheer problem size as- 
sociated with atomistic deconstruction of experimentally rel- 
evant dimensions, typically hundreds of nanometers, as well 
as the scattering physics at the atomistically rough edges. 
We employ a recursive technique to simplify the problem 
size in order to make such a deconstruction tractable. 



2.1 Recursive Green's function Algorithm (RGFA) 

The central entity for quantum kinetics in a weakly interact 
ing system is the retarded Green's function, defined as, 

-1 



' = {EI-H-Ei -E2y 



(1) 



H is the Hamiltonian matrix of graphene, described in this 
paper using a minimal one orbital basis per carbon atom 
with ^0 = — 3eV being the hopping parameter. Ei^2 are the 
self energy matrices from the source and drain contacts, and 
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Fig. 1 Recursive Green's Function Algorithm (RGFA) requires par- 
titioning the structure in left contact, device (layers of N) and right 
contact. 



ri^2 are the corresponding anti-Hermitian parts representing 
the energy level broadening associated with charge injection 
and removal. The corresponding spectral function and quan- 
tum mechanical total transmission are given by 



A = 



(2) 



The local density of states is obtained from the diagonal el- 
ements of A divided by 2;r. 

The contact self energy matrices Z\^2 are calculated from, 
Z = TgT\ where g is the surface Green's function of the 
contact and T is the coupling matrix between the contact 
and graphene. The surface Green's function g is found by 
solving the following equation recursively 



done by stepping through the above equation from right to 
left starting at layer N —I, with the right connected Green's 
function for the Ni\\ layer 

gN^M = [EI-HL-UL-L2]-^ (5) 

Once we reach layer 1 and extract g2,2, the device Green's 
function component for that layer is then calculated from 

^1,1 = {EI-Dl-Zl-h^2g2^2t2^)-' (6) 

where D\ = H\^\^U\. The remaining layer block diagonal 
matrices are calculated by stepping through from left to right 
using, 

^L,L = gL,L + 8L,LtL,L- 1 ^L- 1 ,L- 1 ^L- 1 ,LgL,L (7) 

withL= (2,...,A^). 

2.1.2 From Green's function to conductance and 
non-equilibrium carrier density 

From the diagonal blocks of the layer Green's functions, 
we can calculate the spectral function 

A^L = i{'^m-'^lL) (8) 
The left connected spectral function for layer L is given by. 



^L,L = ^L,l^l,l^L,l 

which requires evaluating a corner block of W 

(10) 

once again by stepping through from left to right, starting 
11 in layer L=2. All the matrices from Eqs. 4pQ 



(9) 



g = [a-TigT{ 
where a = EI 



•ti-i 



(3) 



H for the layered contact unit cell. For 
speedy convergence of g, we use a decimation technique 
outlined in |28|. 

Inverting the matrix in Eq.[T]is computationally expen- 
sive and goes out of hand even for a small structure (width>5nm)!^ 
The Recursive Green's Function Algorithm (RGFA) is a fast 
practical method to avoid inverting the entire matrix by brute 
force, but instead extract only select, relevant blocks of the 
^ matrix, as described in Ref. ||29j, jSOj. For instance the 
transmission across a layered device involves only a corner 
block of the ^ matrix, ^i^l, while the density of states in- 
volves only a diagonal block enabling such partial in- 
version. 



with ' 



are 

NlxNl, where A^^ is the number of atoms in one layer. Since 
matrix inversion is an procedure, inverting substantially 
smaller matrices at a time leads to considerable computa- 
tional savings overall. 

The total conductance is calculated from the transmis- 
sion over all modes 



G = Go-rr(ri,i[Ai,i 



1,1^1,1^1 1 



(11) 



where Go = 2q^/h takes care of spin degeneracy, while val- 
ley degeneracy is captured with the graphene Hamiltonian 

The non-equilibrium carrier density at each layer is cal- 
culated from the above quantities as 
dE 



2.1.1 Evaluating the Green's functions piece by piece 

We consider the device area in layers of N as shown in 
Fig. [T] The right connected Green's function pQ| is calcu- 
lated from 

gL,L = {EI-HL-UL-tL^L^l +gL+l,L+l^L+l,L)"^ (4) 

where Hl is the Hamiltonian of the L th layer, is the 

coupling matrix between adjacent layers and Ul is the elec- 
trostatic potential of the Lth layer. The calculation of gL,L is 



/oo JI7 
^^2^'-[/iAi,L + /2(AL,L-AtJ] (12) 

where /i and /2 are Fermi functions at the source and drain 
respectively. 



2.1.3 Extracting current density 

The current from /th atom to jth atom is calculated from 



Ul 



lij = y I dEIm[^^j{E)Hj^i-Hij<^li{E) 



(13) 
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where the electron correlation function, = '^1}^'^^ and 
scattering function, 1}^ = ^sfs^^Dfu- The source and drain 
Fermi levels are at /i^- = and /i/) = —qVu- To see the current 
distribution in the device, we apply a small drain bias Vd. /^ j- 
is nonzero only if the /th atom and jYh atom are neighbors. 
The total current at one atomic site can be found by adding 
all the components vectorially, // = | Lj/ij'^^^l^ where is 
the angle of the individual current vectors with respect to a 
reference direction. As expected, the NEGF expression for 
satisfies the steady-state Kirchhoff's law, = 0. The 
terminal current Ij is the sum of all currents in one atomic 
layer, and should be equal to the current from the Landauer 
formalism using the total transmission. 



, q 



[fs{E)-fo{E)]dE 



(14) 



While the total current It (Eq. [14]) can be calculated effi- 
ciently using the partial blocks in the RGFA, the spatially 



resolved current density (Eq. 13), useful primarily for vi- 
sualization purposes, employs the full matrix ^ and is thus 
computationally expensive. The recursive algorithm to cal- 
culate blocks of is described in |29j . The current density 
from layer L to L + 1 is given by, 

h^L+l{E) = ^Im[tL,L+l'^E,L+liE)-tL+l,L'^E+l,Lm(i5) 

The matrix (E) has the size of A^^ x Nl carrying the 

atom to atom current components between the two layers. 

2.2 k space formalism (KSF) 

We can further cut down the problem size when the sys- 
tem is periodic along the direction perpendicular to electron 
transport, by employing B loch's theorem. The periodicity 
allows us to decompose the system so that the transverse k- 
points act as decoupled 1-D chains or modes, whose contri- 
bution to the total current can be calculated independently. 
We can thereafter calculate the conductance through bulk 
graphene or graphene nanoribbons (GNRs) with mode by 
mode resolution. As we shall shortly see, it also allows us to 
visualize the transition from the smooth transmission pro- 
file of bulk graphene to the stepwise transmission of GNRs 
driven by width quantization. 

Consider a graphene sheet with an applied potential vary- 
ing only in the x direction (Fig.|2^). After dividing the sheet 
into blocks, we can split the hamiltonian into 4x4 block 

matrices (H p g) describing the interaction between the wave 

yi yj 

functions x^f p and i/A a belonging to blocks and j^, with 
subscripts denoting the rows labeled along the transverse 
direction, and superscripts denoting layer columns labeled 
along the transport direction (Fig. [2]b). Using this pictorial 
representation of the Hamiltonian, the Schrodinger equation 
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Fig. 2 a) Graphene sheet split into blocks that reflect its periodicity 
in the y direction, b) Representation of Graphene 's hamiltonian split 
into matrices describing the interaction between block's wave func- 
tions. Each box represents the interaction of the a block with itself 
and each line represents the interaction of a block with a neighboring 
block, c) After Fourier transforming the Graphene sheet can be split 
into independent chains, one for each ky. 



for the yf block is given by 

E Xj/v =Hv p xj/ p + Hpv wp + h\ p xj/ p 

^ yi yi yi ^ yi yi yi+i ^ yi+i }^f_i}'f 

-^^yfyr'i'yr' ^Kr'y^'^y^^'' ^^^^ 

To take advantage of the manifest periodicity in _y, we as- 

ik ^ 

sume Bloch type solutions x^f^q = xj/^e ^'^j , with y'j the vec- 
tor describing the position of the block y^j. Substituting this 
solution into Eq. [16) we get 



with 



II 



P,P+i 



- H p p+i , 

ytyi 



^ikyb ^ 



H 



yf-iyi 



TtP-^,P _ TT 

My — rt p-i p. 

h yi y'i 



(18) 



(19) 



(20) 



For each ky, Eq. 17 with a varying index p can be interpreted 
as a one dimensional chain (Fig. [2]:) decoupled from chains 
with other ^^s. In other words, we have split the graphene 
sheet into a set of independent chains or modes. Now, we 
can use the usual NEGF formalism to find the transmission 
of each chain or mode (Ti). 



ny = Tr{r,%r,l<=^l) 



(21) 



The quantities (E, F, G and T) in real space are found by 
inverse Fourier transform, since our Bloch assumption is 
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equivalent to Fourier transforming from the discrete space 
to Z^^;- space. For example, the Green's function is given by 

^(3f.3'P = ^I%^"^^^'^"'^^ (22) 

with Nk, the number of points in ky space, and the ky depen- 
dent Green's function given by 

^ky = [{E^iri)I-Hk^-E,^ky-^2^ky]-' (23) 

The sigma matrices are calculated as before, E^^ = tJ^^^Ti 
but the surface Green's function is now ky dependent. 

gky = [OCky-^lgky^lr^ (24) 

where a^^ = EI — H^^ . H^^ takes care of all the couplings in 
the transverse direction. 

Thus the combination of semi-infinite contacts in the ±x 
direction and Fourier transformation along the ±_y direction 
amounts to a change in basis set for the completely periodic 
bulk 2D graphene. Since trace is preserved under the basis 
transformation, we can calculate transmission directly in ky 
space 

T = Y,Tky DOS=^Tr{^-^^) (25) 

ky 

where 

'^ = Y.% (26) 

k 

The distinction between bulk graphene and GNR arises pri- 
marily from the different sets ofky^ (in practice, one needs 
to worry about edge reconstruction, next nearest neighbor 
interactions, strain and roughness at the edges of GNRs. For 
an ideal GNR of width W =Nb= N^b the hard wall bound- 
ary conditions restrict ky to the discrete set 

nTt 

ky = ^ with < kyb < 71 and n integer. 

In the case of bulk graphene we can assume periodic bound- 
ary conditions every W = Nb, which restrict the ^3;'s to 

2n7T . , ^ ^ , . 

ky = with — TT < kyb < 71 and n integer, 

and let A/^ ^ oo. Note that as N increases, the transmission 
evolves from a stepwise curve showing width quantization 
to a smooth curve. 

2.3 Applying open boundary conditions 

Since simulations deal with finite sized matrices, it is 
important to make sure the results do not get influenced by 
edge effects. While realistic physical devices do have edges. 



their impact on transport is complicated, dominated by strain 
and reconstruction, roughness, possible charge trapping, edge 
dipoles, localized vibrons and other subtleties. In fact, mea- 
surements of tilt-dependent junction resistance in graphene 
show a surprisingly weak influence of edge scattering 1 16j 
[STI . For modeling purposes therefore, it is important to en- 
sure that additional edge effects do not creep in simply be- 
cause of the finite- sized domain of our simulation grid. 

Section 2.1 deals with hard- wall boundary conditions 
along the transverse direction. For such structures, elimi- 
nating edge effects requires making the widths very long 
(often longer than device length), which may have added 
repercussions on material properties (see section 3.1). Sec- 
tion 2.2 implemented periodic boundary conditions (PBC) 
instead with a k-space formalism, allowing us to work with 
just a couple of nearest neighbor unit cells along the trans- 
verse direction, but requiring strict periodicity along that di- 
rection. In this section, we will discuss how to implement 
open boundary conditions (OBC) that would allow us to ir- 
reversibly remove an electron impingent at the edges. The 
trick is to treat the edges as virtual 'contacts' with self-energy 
matrices built out of complex numbers, acting like energy- 
dependent lossy potentials. 

In this section, we will show how to generate exact open 
boundary conditions (OBC) at the edges of a device. It is 
worth emphasizing at this point that these are quantum bound- 
ary conditions, i.e., acting on the retarded Green's function 
^. What is much harder to solve is a thermal open bound- 
ary condition acting on because this depends on details 
of the scattering processes outside the simulation regime. 
In metallurgical contacts, we impose by fiat Fermi-Dirac 
distributions, but in the extended geometry just below our 
simulation regime, the distribution is non-equilibrium and 
the corresponding boundary conditions need to be estimated 
self-consistently with the inside the device in order to 
prevent the outside regime from sucking out the electrons 
too aggressively or else inadequately. 

2.3.1 The 'scooped' open boundary self -energy 

The 'obvious' way to implement some kind of open bound- 
ary is to attach virtual leads along the transverse direction 
and extract their self-energies ^3^4 using recursion. This is 
schematically described by Fig.[3]^h). While this method will 
take the electrons away, the leads do not completely span the 
outside regime relative to the central device represented by 
the Hamiltonian [H], and the missing chunks at the corners 
are expected to create spurious reflections. One could extend 
the contacts L\^2 or ^3^4 laterally using k space formalism to 
span those comer regions, but the problem is that the four 
virtual leads are ultimately decoupled from each other with 
no bonds running in between, and thus do not quite represent 
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(9) (h) 

Fig. 3 Tracking the injected Gaussian wavepacket at different times (a-e), the wavepacket escapes due to the exact open boundary condition (g) 
imposed through the self energy term. We see some reflections along the corner if we adopt a less exact (but computationally less expensive) 
method to impose open boundary condition through transverse contacts (h). 



a straightforward extension of the device domain beyond the 
simulation regime. 



What a proper OBC will need is a geometry like Fig.[3jg), 
where the virtual 'contact' should be a single monolithic 
structure with the central device region 'scooped' out. The 
single self-energy matrix Eopen for this outside regime must 
be reverse engineered so that the Green's function of the 
central region matches that of the infinite system periodi- 
cally extended along all axes. Once we find this self energy 
(section 2.3.2), we can use it for any modified condition (e.g. 
a gate voltage, a molecular adsorbate or an injecting contact) 
that only influences the central region but keeps the outside 



unaltered. 



2,3,2 Calculating self energy E^ 



'open 



To find the boundary free real-space Green's function in 
the central region, we first extract the k space version corre- 
sponding to the infinite periodic system. 



= [(^ + iri )I-Hy,- E, ,k - i:2,k] " 



(27) 



where the k-dependent quantities are Fourier transformed 
sums over nearest neighbors. From this we can calculate the 
real-space Green's function projected onto the device region 
by inverse transforming 



:£^kexpHk-(R„ 
k 



■Rn)]/N,nNn 



(28) 



for the atoms 'm' and 'n' spanning the central device seg- 
ment. For a system with A' atoms, is N x N. 
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We can now reverse engineer the open boundary self- 
energy, armed with this real space Green's function and the 
N xN device Hamiltonian H, giving us 



(29) 



This self energy can be used to calculate transmission through 
the structure with open boundary conditions, additional in- 
jecting/removing leads El ^2 and any other perturbations act- 
ing only on the central region through a localized potential 
V that does not influence the boundaries 



Trace (Ts^new^D'^lew ) 

{EI-H-V-Ls-Ld 



^open ) 



(30) 



2,3.3 Verifying the open boundary condition 

In Fig. [3] we verify the efficiency of the open boundary 
condition by launching a Gaussian wavevector with a spread 
{(7x, CTy} and initial quasi-momentum set by the wave- vector 

^'{x,y,t = 0)= ^ 



27r y^ojoy 



(31) 



about an injection point (xo,jo)- We calculate the evolu- 
tion of the gaussian wave-packet with time by numerically 
solving the time-dependent Schrodinger equation within the 
Crank Nicholson algorithm (described below), acting on the 
time-dependent vector {^{t)} with row entries correspond- 
ing to the spatial coordinates 

{W{t + At)}=[M]{W{t)} 

1 -1 



M 



'open 



At/2 



I^i{H^Zopen)At/2 



(32) 

In the absence of any open boundaries, the wavepacket stays 
confined within the simulation domain and bounces around 
(not shown). With the scooped self-energy I^open (Fig. [5), 
the open boundary condition takes the electron out at the 
boundary (Figs. |3^-e). Note that the self-energy we used 
works only at a specific energy E set by the pair (kj^ky), 
although we can extend it to all energies by Fourier trans- 
forming. 

Fig. |3jf) shows the results for a cross-geometry with 
transverse contacts represented by Fig. [3jh). As expected, 
the missing pieces along the corners give a sizeable reflec- 
tion of the wavepacket. 

3 NEGF simulation of electron transport in graphene 

In the following sections, we will use a combination of nu- 
merical techniques outlined so far to simulate current flow 



through graphene segments comparable in size to experi- 
mental dimensions. We show how these simulations accu- 
rately capture the nuances of conductivity, conductance, mag- 
netotransport, pseudospintronics and Klein tunneling within 
a common, unified simulation platform. 

3.1 Minimum conductivity vs conductance 

One of the unique properties of graphene is that its mini- 
mum conductivity G (rather than conductance G = gW/L) 
is quantized in units of Aq^ /nh, when the width to length 
aspect ratio W/L is large. In other words, G is proportion- 
tal to the aspect ratio W /L 'm wide graphene samples. The 
width-dependence is expected in a large ballistic device as 
the number of modes is roughly proportional to the number 
of Fermi half- wavelengths one can fit into W . However, the 
average transmission per mode, and thus the conductance, 
is usually length-independent. This is because propagating 
modes have transmissions of unity in a ballistic channel, 
while evanescent tunneling modes have transmissions that 
are nearly zero. However since graphene is semi-metallic, 
a wide sheet supports a nearly continuum set of sub-bands 
with ultralow band-gaps and tunneling probabilities that are 
not ignorable. In fact the ~ 1 / cosh^ qyL tunnel transmis- 
sion terms {qy = nn/W) from these closely spaced modes 
all add up to an overall 1/L dependence of the total trans- 
mission (more precisely, a factor 2W/L71) and therefore the 
conductance G, leading to the conductivity quantization. Ex- 
periments on dirty samples show such a quantized mini- 
mum conductivity, albeit off by a factor of ~ 3 — 4, usu- 
ally attributed to Coulomb scattering from charge puddles 

We will use the NEGF technique to extract the conduc- 
tivity of graphene at the Dirac point. To recap, for large scale 
graphene devices (L > 1/im), the conductance per mode is 
expected to be Aq^ /h (the mode count is given by the inte- 
ger part of Wkp / tt), while for smaller lengths (L < 500^m) 
and larger widths, the conductivity is expected to approach 
Aq^ / {nh). We will simulate the structure shown in Fig.jJJa) 
to calculate the minimum conductivity of a clean graphene 
sheet. The graphene device is connected to two heavily doped 
graphene contacts. As long as the contacts provide a large 
number of modes and thus minimal series resistance, the de- 
tails of the contact metallicity should not matter. The en- 
ergy band diagram for the problem is shown in Fig. ^b). 
We model the system atomistically with RGFA as described 
in section 2.1. 



We use Eq.[TT]to calculate the total conductance G. Min- 
imum conductivity Gmin = GL/W is shown in Fig.[4jc) along 
with the experimental data from Ref. p3| . We take two dif- 
ferent widths W = 100, 200 nm to calculate G and Gmin for 
various lengths (L). Gmin is found to be very close to Aq^ / 7th. 
As expected, the minimum conductivity stays close to this 
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Fig. 5 (a) Density of states (DOS) and conductance (inset, in units 
of 2q^ /h) calculation from atomistic tight binding RGFA for two 
different widths, (b) Similar calculations over a wider energy range 
done with KSF, atomistic NEGF within RGFA, and an analytical in- 
tegration of E — k dispersion. Conductance (inset) is in units of 
2q^ /h X 1000/im~^ Conductance from KSF (red line) agrees with sim- 
ple linear approximation of no. of modes (black circles) at low energy. 



Fig. 4 Numerical extraction of graphene minimum conductivity in the 
small dimension limit (L<500nm). (a) Schematic of the device, (b) 
Energy band diagram for the device with doped contacts, (c) Calculated 
Omin and (d) G^/^for various aspect ratios along with experimental data 
fromRef. |33|. 

number as long as W /L is large and a continuum evalua- 
tion of the tunnel contribution is a valid approach. In the 
opposite limit W <C L the conductivity increases with length 
and we reach the conventional regime of conductance quan- 
tization (Fig. |4jd)). Note that for a graphene ribbon of fi- 
nite width, the valley degeneracy of the fundamental mode 
is broken f3?|, so that the minimum conductance is 2q^ /h. 
At higher gate voltages as more modes come into play, the 
conductance per mode reaches 4q^/h, with the degeneracy 
of 4 coming from the presence of 2 spins and 2 valleys at the 
same energy. 



3.2 RGFA for finite sized graphene devices 

We will now apply the RGFA technique discussed in 
section II to simulate three types of finite sized graphene 
based devices - uniformly doped, pn step junction and npn 
barrier junction. 

3.2.1 Uniformly gated graphene 

We apply RGFA to calculate the density of states (DOS) 
and conductance of graphene sheets of various widths (W). 
Fig.jSja) shows the DOS (and conductance in inset) for W = 
10 and lOOnm. Narrower nanoribbons show a gap as well as 
Van Hove singularities due to quantization for select chiral- 
ities while wider sheets show a quasi-linear bulk graphene 
density of states. In practice, nanoribbons also need atten- 
tion to edge state dynamics, particularly the presence of strain 
and roughness |T2| . 

A striking property of graphene is its anomalous inte- 
ger quantum Hall effect. When the Fermi energy lies be- 



tween two Landau levels (LL) in presence of a magnetic 
field, conventional two-dimensional electron gases show a 
vanishing longitudinal resistance p^x = and a quantized 
Hall conductance Gh = 2q^N /h where A/^ is a non-negative 
integer representing the number of filled Landau levels. In 
contrast, chiral quasiparticles in graphene exhibiting Berry 
phase show | 36p7 | a Hall conductance G// = ^ (A^+ 1/2), 

^ ^inr.(r.\ . whcrC Tc = • 



with LL defined as En = ^sign{n) 
is the cyclotron radius. This non-conventional sequence can 
be explained with the presence of a LL at E = resulting 
in a four fold degeneracy at zero carrier density (from spin 
and valley degeneracy). On the contrary, Bilayer graphene 
LLs are defined as, — h(Dc\/n{n — 1), where (Dc is the 
cyclotron frequency. Now both E^ and E\ are at zero energy 
producing thereby an overall eightfold degeneracy at zero 
carrier density and QHE plateaus of G// ^(A^+1) jSSf. 

Numerical modeling magnetotransport in graphene re- 
quires a minor modification to the transport scheme outlined 
earlier. We replace the kinematic momentum with the quasi- 
momentum k ^ k — qKjh, so that the plane wave terms 
in the Bloch representation pick up an additional phase of 
^-^^/^/A di ^ magnetic vector potential such 

that V X A = B. Thus the hopping parameters between atoms 
'm' and 'n' are now given by 



:^oexp 



rn 
Jm 



A dl 



(33) 



For z directed magnetic field, our guage is, A = {—By ^0^0). 
We can now turn on a magnetic field perpendicular to the 
sheet, modify the hopping terms as above, and extract the 
transvese (Hall) conductance as a function of varying Fermi 
energy location. The Hall conductance Gh = I/Vh, where 
/ = {2q/h)N{ilL — {Ir) is the current, N is the number of 
filled Landau levels and the Hall voltage Vh = qx the differ- 
ence between the electrochemical potentials of the -\-k and 
—k edge states. These states represent skipping orbits along 
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Fig. 6 (a) In presence of a magnetic field, the current is carried by the 
edge states that separate into -\-k states at the upper edge in equihbrium 
with the left contact and —k states at the lower edge in equilibrium 
with the right contact, (b) The Hall voltage goes through plateaus at 
^ (A/^+ 1/2) for single layer and ^ (A/^+ 1) for bilayer graphene with 
N a non-negative integer. Note the presence of a jump at zero energy, 
which is not seen in a two-dimensional free electron gas (and arises 
from a half-filled Landau level at the Dirac point), and the additional 
factor of 2 in the plateau heights, arising from valley degeneracy, (c) 
For a placement of the Fermi energy between two Landau levels, elec- 
tron distribution function at the top and bottom edges resemble Fermi- 
Dirac distributions of the corresponding contacts. Since each current 
carrying state sees a constant electrochemical potential (along trans- 
port direction x), the longitudinal resistance vanishes. 



the edges created by the cyclotron orbits in the bulk, and 
are separately in equilibrium with the left and right contact 
Fermi energies. 

Fig.jSJ) shows the calculated Hall conductance yielding a 
series of plateaus for both single layer and bilayer graphene. 
Fig. [6]: shows the local carrier distribution functions /(±k) 
obtained from the ratio of the local carrier density and the 
local density of states, in other words, the ratio 



/(^)=^L,L(/,OAm,L(/,/)-^L,L(/,/)) 



(34) 



where / is the index for an atom belonging to either top or 
bottom edge. The Hall voltage turns out to be V// = ^(/i^ — 
IIr) (Fig. [6]:) giving conductance plateaus for Gh- The elec- 
trochemical potentials (iiljI^r) don't change along transport 
direction x so that the longitudinal resistance (given by the 



drop in electrochemical potential along the channel) is zero. 

3,2,2 Single p-n junction 

junction heterostructures in graphene act qualitatively 
different from conventional pn junctions. In normal semi- 
conductors, the presence of a band-gap blocks electron flow, 
so that reducing (increasing) the built-in potential across a 
pn junction with a drain bias leads to an exponentially in- 
creased (reduced) current, creating a rectifying current- voltage 
(I-V) characteristic. The lack of a bandgap removes such 
rectification in graphene pn junctions (GPNJs). However, 
the chiral nature of the electron flow provides fascinating 
electron 'optics' behavior, reminiscent of Snell's law, albeit 
with some notable twists. Across a junction for instance, 
the conservation of transverse quasi-momentum {ky) leads 
to anomalous Snell's law 



Kp\sinB\ = —Kp2sinQ2 



(35) 



where the refractive index is set by the Fermi wavevector 
and thus the gate voltage applied. The opposite sign of the 
voltage across a pn junction, arising when the quasimomen- 
tum component flips sign in going from conduction to va- 
lence band, means that the system acts like a negative index 
metamaterial. 

While the electron trajectories are set by the above Snell's 
law, i.e., the conservation of transverse quasi-momentum, 
the actual transmission probabilities are set by the conserva- 
tion of pseudospins. The pseudospins arise from the eigen- 
vectors corresponding to the photon-like graphene E-k. The 
eigenfunctions are given by 



1 



se 



'6 



(36) 



where is the 2-D angle setting the quasi-momentum k and 
5" = zbl denotes the conduction or valence band (An addi- 
tional sign flip with respect to kx occurs near the other val- 
ley). Up and down pseudospin states are given by = 
and 7t, corresponding to the bonding and antibonding com- 
binations of the Pz dimer basis sets in graphene. Analogous 
to optics, we can derive the equivalent 'Fresnel equation' 
for transmission probabilities, by matching the pseudospinor 
components across the junction. Assuming a split distance 
'd' between the backgates defining the pn junction, and a 
voltage difference Vq between the gates, the transmission 
works out to be 



cos 6i cos 62 
61^62 



-nhvFkpdsin^ Gi /Vq ^^j^ 



below a critical angle defined from the Snell's law (Eq.[35j, 
61 (62) representing the incident (refracted) angle for elec- 
trons. The first term shows that at normal incidence (61 =0) 
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Fig. 7 Numerical calculation (NEGF) of a single pn junction conductance for both abrupt and smooth GPNJ. (a-b) Variation of total resistance 
with dopings in the two regions for abrupt (a) and smooth junction (c) Resistance as a function of Fermi energy {Ep) for a fixed built in potential 
Vb = 0.2eV. The red (abrupt) and black (smooth junction) lines are resistance plots along the corresponding black and red lines in (a) and (b). 
Conductance (inset) which can be viewed as G = G^MTav where Tav, pinches off at two points - one due to vanishing M and at another one 
due to vanishing Tav, average transmission per mode, (d) The resistance asymmetry between pn and nn regime from plots along specific doping 
{AE\ lines from (b). (e) Extraction of Tav numerically. The transmission of a symmetric GPNJ is 2/3 (at Ep = 0), which results in asymmetric 
conductance vs doping in GPNJ. The solid lines are from NEGF and the circles are from analytical calculation, Eq.|^ (f) Junction resistance (Eq. 
[40) enhancement due to a tilt in the junction with respect to the tilt. This time we vary the doping along the diagonal (orange line in (b)) as done 
in the experiment |31 1. 



the transmission function is unity regardless of the size of 
the barrier Vq. The absence of normal reflection is a manifes- 
tation of Klein tunneling, and arises because backscattering 
requires the flipping of pseudospins, in other words, a fast 
spatially varying potential. The second, exponential factor 
(drops out for unipolar n^n or p^p junctions) in the trans- 
mission equation comes from tunneling across the slowly 
varying junction, where the transverse modes face a momentum- 
dependent band-gap, reminiscent of cut-off frequencies in a 
wave-guide. 

The overall conductance of a GPNJ can be calculated 
by summing over all transmissions, G = Gq YJq ^i^)' which 
can be thought of as the total number of modes times an 
average transmission. Thus the average transmission over 
all modes at a particular energy is given by Tav = G/GqM, 
where number of modes, M = A\Ep\W /nfivp. From RGFA 
calculation, we get the total conductance G = rr(ri^J2^^) 
in units of Go = 2q^/h. 

Figs. |7] (a,b) show the doping-dependent resistance of 
a GPNJ for a lOOnm wide graphene sheet, for abrupt and 
slowly varying junctions respectively. In fig. [7] we plot the 
variables against the shifts in the Dirac point in both regions, 
i.e. A El and AE2. The resistances at the upper left and lower 
right comers of the plot are higher than the other two (pn vs 
uniformly doped), resulting in an asymmetric resistance vs 
doping in fig. |7jd) (plotted for specific doping values AEi 
as indicated with horizontal lines in|7jb)). 



The WKB term in Eq.[37]is only present in the pn junc- 
tion regime, and that is why only the pn junction resistance 
is affected while going from abrupt to smooth junctions (a- 
b). For di pp or nn junction the Fermi energy does not cross 
the smoothly varying Dirac point anywhere in the device and 
the transmission expression only includes the wave-function 
mismatch term. Fig. [tJc) shows the resistance variation for 
a fixed built-in potential Vq^ AE\— AE2 (along black and 
red lines in fig. |7ja) and (b)). In the conductance variation 
(inset), we see effectively two Dirac points, which has a sim- 
ple physical explanation. Recall that the normalized conduc- 
tance can be decoupled into the mode count from one end 
times the average transmission over to the other side 

G/Go=MiTn=M2T2i (38) 

The left conductance minimum at Ef = —O.leV is where 
Ml becomes zero, while the right one at Ep = -\-0.1eV is 
where the average transmission for all modes, 7i2, becomes 
zero. To make this clear, we calculate Tav numerically (Fig. 
[TJe)), clearly showing a vanishing transmission at the second 
Dirac point. To calculate Tav, we first simulate a graphene 
device with uniform doping and extract the overall mode 
count from the ballistic conductance (Gi = MGq). We then 
simulate the device with different dopings (finite built-in po- 
tential, Vb, and conductance G2 = GoMTav)- The ratio of 
G2{Ep) and G\{Ef) at each energy gives us Tav{Ep). For 
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X(nm) X(nin) 

Fig. 8 Tracking electron trajectories in GPNJ, negative index makes the refraction angle negative - resulting in spatial focusing of current flow 
for equal dopings (left column). Right column shows total internal reflection when the refracted side has lower doping, thus reflecting electrons 
outside the critical angle. 



example, for a symmetric GPNJ (at Ep 
T{e) 



0), we have 



/TT/Z 

r7t/2 
J-71/2 



AO 



-de 



cos^e 



kpcosGdO — 



WkF /"^/^ 



271 



rK/Z 

/ ^ 

J-k/2 



cos^BdB 



-M(E) 



(39) 



leading to T^y = 2/3, which is close to what we get in the 
numerical calculation. Thus the conductance of a symmet- 
ric GPNJ is 2/3 of the uniformly doped graphene sheet with 
the same Fermi energy (Ep). In Fig.|7je) we also show the 
analytical calculation of Tav (circles) in the same manner as 



done above but this time using the transmission in Eq. 37 



The analytical calculation matches very closely with the nu- 
merical calculation (solid lines). 

A more direct evidence of chiral tunneling is seen from 
the resistance across a tilted junction. We use the follow- 
ing equation to extract the junction resistance |7 |, that elim- 
inates the quantized resistance from the contacts. 



Rj 



i'4r 



■ MTav 



(40) 



We calculate M and Tav numerically from NEGF as dis- 
cribed earlier. Fig.[7]f shows the variation of Rj plotted against 
AE\ (along orange line) for a lOOnm wide graphene sheet. 
For low carrier density in the middle, the background dop- 
ing (^o) makes the device an ri^n or p^p junction and 7?/ is 
low. For n\^2 > ^o, i-e., in the wings of the plot, the device 
reaches the — n junction regime and R j becomes high. This 
characteristic doping dependence of junction resistance has 



been seen in experiments in the past |[39j|40|, usually cap- 
tured by the odd resistance, 2Rodd = R{^i ^ ^i) — R{^\ ^ —^i)- 
The chiral tunneling can be seen more directly by varying 
the tilt angle, as seen in experiments recently 1 16 |. The pur- 
ple squares show the trend for a tilted GPNJ, where a pro- 
nounced peak is seen mRj. With increasing tilt between the 
junction and the transport axis normal to the contacts, the 
incident angless of the incoming modes are effectively in- 
creased, leading to higher resistance. Such an increase in 
junction resistance is seen in a recent experiment | 31 1 and a 
detailed theoretical treatement can be found in 1 16 |. Interest- 
ingly, the results suggest the absence of specular reflection 
from the edges of the graphene sample. 

We will next apply the current density formalism in RGFA 
to show the electron trajectories inside the GPNJ device. 
The total current at an atomic site is equal to vector sum 
of all current components to nearest neighbor atoms (from 



Eq. 15 ). The classical ray tracing analysis from Snell's law is 
shown in Fig.[8ja). For equal dopings on both sides {kF2 = 
kpi), the angle of refraction is exactly equal (with a neg- 
ative sign) to the incident angle. As a result, a group of 
electrons originating from a point contact focus back to one 
point on the refracted side |6 |. Fig. [8jb) shows the trajec- 
tories when the doping at the refracted side is smaller than 
that at the incident side {kF2 < ^fi) making the critical an- 
gle smaller than 7r/2. The geometrical 'optics' trajectories 
corresponding to electron focusing and total internal reflec- 
tion are reproduced by the calculated current density with- 
ing an atomistic RGFA simulation, shown in Fig. [Sjc), (d). 
A small source contact (lOnm wide) is placed 50nm to the 
left of the pn junction, and electrons are injected with a 
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small drain bias (Yd = 0.08 V) around a Fermi energy. When 
the gates are biased symmetrically around the junction, al- 
though the electrons see a voltage bias along the drain that 
spans the entire device width, the pn junction Hamiltonian 
and associated pseudospin conservation causes the electrons 
to focus to a small point at the drain. For an n+p junction, 
the electrons incident above the critical angle are unable to 
preserve their transverse quasi-momentum and reflect back, 
while those within critical angle tunnel to the opposite band 
on the other side. 

We now investigate the quantum Hall plateaus for a sin- 
gle pn junction. According to recent experiment |41 1, the 
conductance in the pn junction regime follows the Ohmic 
conductance rule 



G{q'/h) : 



Ohmic plateaus 



(41) 



Vyi and Vp are the filling factors in the two segments of the 
device. This was explained by Abanin et. al. | [42} in terms 
of mode mixing at the p — n interface for large systems with 
diffusive transport. As graphene filling factors are {2, 6, 10, ...}, 
it produces junction plateaus as {l, |, ...^^c.} for filling 
factor combinations of (2,2), (2,6) etc. matching closely 
with the experiment | 41 1. 

This result changes considerably for coherent ballistic 
transport for smaller structures. Tworzydlo et. al | 43 1 showed 
that for armchair nanoribbons the conductance plateaus be- 
come 

G = ^Go(l — cos<P)\ Ballistic plateaus (42) 

independent of individual filling factors v„, v^. Go = 2q^ /h 
is the lowest Hall plateau for graphene. <P is the angle be- 
tween valley isospins at the top and bottom edges and de- 
pends on the chirality of the nanoribbon. Depending on the 
number of hexagons N along transverse direction, = tt for 
N = 3M+ 1 and <P = 7l/3 otherwise. This leads to 



//N = 3M+1 
otherwise 



(43) 



The ballistic to ohmic cross over for the plateaus can be re- 
covered by using interface and edge disorder |44|. 

The magnetotransport for zigzag ribbons cannot be ex- 
plained with valley isospin arguments. Akhmerov et. al. | [45l 
proposed the theory of valley valve effect with intervalley 
scattering which leads to complete suppression of conduc- 
tance in the zigzag /injunction, if the number of atoms across 
the ribbon is even. 



G{qVh) 



if N = 2M 
otherwise 



(44) 



Impact of depletion width d (spacing between two gates) 
on GPNJ ballistic magneto-conductance for both armchair 
and zigzag ribbons are shown in Fig. |9] Magnetic field B 




10 20 30 40 
Depletion width (nm) 



10 20 30 40 50 
Depletion width (nm) 



Fig. 9 Magnetotransport in graphene pn junction. Magnetic field, 
B=10T and filling factors are Vn = yp = 6. (Left) Change of conduc- 
tance with depletion width, d for Armchair ribbons. For high depletion 
width, we get the expected plateaus of 2 and ^ . (Right) Plateaus of 
and 2 for Zigzag ribbons. Dotted horizontal lines show analytical pre- 
dictions. 




Fig. 10 Physics of npn junction, (a) Band diagram. The black line 
shows the change of Dirac point as a result of different dopings at dif- 
ferent portions of the device. Blue (red) region in the E-k indicates 
empty (filled) states, (b) Electron trajectories in such junctions, mul- 
tiple focusing takes place, (c) Ballistic conductance as a function of 
Fermi energy for a fixed barrier height for different barrier width, D=5, 
25 and lOOnm (black, red and blue) . Width of the device is W = 50nm. 
The solid lines are from numerical NEGF calculation, while the cir- 
cles are from analytical calculation, Eq.[45] (d) Experimental data from 
Ref. | 46 1. Conductance in the experiment is much lower due to scatter- 
ing mechanisms. 



= lOT and dopings are AEi = AE2 = 0.1 5eV producing 
= Vp = 6. Conductance values approach analytically pre- 
dicted numbers only when the depletion width is sufficiently 
large (> 30nm), filtering out higher Landau Levels (LL) 
| [44| . It requires larger depletion width as individual filling 
factors get higher. In most cases, the lower plateaus (^ and 
for armchair and zigzag respectively) are more difficult to 
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achieve than the higher one (2). 



Single layer 



Bilayer 



3,2,3 n-p-n junction 

By applying a step Hke potential, we can realize an npn 
junction. The dopings can be done by either a combination 
of global back gate and top gate | 40 | or by selective chemi- 



cal doping. In fig.[TO[a), we show the energy band diagram 
for the device. We consider both junctions as abrupt, J = 0, 
while the extent of the p type barrier region is D. In fig. 
7Q|c), we show an NEGF calculation of the npn conductance 
for various D values. For larger D, two Dirac points are ev- 
ident from the pinched off conductance, but for smaller Ds 
we see a considerable tunneling around the second Dirac 
point (Ef = 0.5eV), where the critical angle for the incident 
electrons is supposed to be very small. The first dip happens 
at the Dirac point of the incident n region where Mi = as 
before, while the second corresponds to the electrons in the 
n region aligned with the Dirac point of the barrier p region 
and the transmission 7i2 is small. In a single p — n junc- 
tion, the modes with higher angles than the critical angle 
reflect back as they do not have any propagating states to 
tunnel into while preserving their transverse quasimomen- 
tum. But in the npn junction case, the length of the forbid- 
den region D is finite, and the electrons even while aligned 
with the Dirac point of the central p region can still preserve 
their quasi-momenta by tunneling to the other side. This in- 
creased conductance has been seen experimentally | [46| in 
the past (Fig. [TO]). In the limit when D is very large, we 
approach the single p — n junction case we worked out in 
Fig. [tJc) inset. And when D is very small, we approach the 
uniform doping (n doped here) case through significant tun- 
neling. Note that we also see oscillations for junctions, 
which originate from the Fabry Perot cavity formed by the 
electrostatic barrier |'46',^7 |. Both the tunneling and the res- 
onant oscillations can be captured analytically by matching 
eigenvector components across the pnp junction, as worked 
out in Ref. iSj. 



2ie^^ sin{qxD) x (sinij) —ss^sinO) 



ss'[e-'^-^cos{(l) + e) + e^i-^cos{^ - B)] 



- 2isin{qxD) 
(45) 



giving us a transmission, r(0) = l—R= 1 — |rp. ^^^^ is the 



wave vector inside the barrier, qx= y {Ef — Vq)^/ {h^vj^) — kj, 

6 = tan~^ (ky/qx) and s = signEF and s' = sign{EF — Vq). 
We then sum over different modes, set by the width of the 
graphene sheet, to get the total conductance G{Ef) in units 
ofGo = 2q^/h (Fig. 
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(c), ,circles), in excellent agreement 
with the NEGF atomistic calculation and qualitatively with 





Fig. 11 Angle resolved transmission in single layer and bilayer 
graphene for a fixed Fermi energy and for various built-in potential 
Vq. At left (Ef = 0.5eV), Klein tunneling with maximum transmission 
at normal incidence. At right (Ep = 0.05eV), Klein reflection in bi- 
layer graphene for normal incidence. Top pannel shows the forward 
and backward moving branches colored according to the pseudospin 
texture for normal incidence. 



experiments (Fig.[TO])). 

3.3 KSF to simulate bulk graphene and GPNJ 

The inset in Fig.[5jb) shows the calculated conductance 
of bulk graphene from KSF. Note that the peaks due to finite 
width quantization in Fig. |5] are removed by the application 
of periodic boundary conditions. The black circles show the 
trend that is expected from a linear approximation of the 
number of modes M = 2W\EF\/7lhvF at low energy. The 
main figure in Fig. [5jb) shows the bulk graphene density 
of states calculated from three methods. The red line shows 
the result from an atomistic NEGF calculation. The black 
circles are from an analytical formula |48|, while the blue 
line is from a numerical extraction of the density of states 
(DOS) from the graphene E — K relation. All three calcula- 
tions match very well with one another. 



KSF approach also allows us to show the angle (mode) 

to 



dependent chiral tunneling in graphene. We use Eq. 21 



get which is converted into a polar plot (Fig. \\A\ using 

e 



- sin~ 



^ ^- for transmission per spin per valley. For single 



layer graphene the transmission becomes unity for normal 
incidence (known as Klein tunneling). On the contrary, for 
bilayer graphene, the reflection becomes unity (Klein reflec- 
tion). The tranmission ultimately depends on how well the 
wave-functions match on both sides of the junction (branches 
with same pseudospin components have same color in Fig. 
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— d = 1uin 
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\ —abrupt GPNJ 
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Fig. 12 KSF-RGFA simulation of smooth GPNJ. (a) Conductance 
characteristics as a function of Fermi energy for a fixed built-in poten- 
tial Vq = 0.6eV, conductance is suppressed in the pn junction regime 
due to the exponentially decaying transmission of higher order modes. 
Inset: Average transmission per mode Tav for a symmetric GPNJ (Ep = 
point in (a)) vs the separation between two gates. 



[TT]). 

3.4 Combining KSF and RGFA 

We now combine the two transport formalisms KSF and 
RGFA to a graphene p — ^junction and show how this method 
allows us to simulate very long, micron-sized devices with 
periodic boundary conditions along the transverse direction. 
We apply this technique to numerically reproduce the tun- 
nel transmission across a GPNJ with split gates, which re- 
quires simulating a long device to extract the proper de- 
cay constant. We assume that the Dirac point varies linearly 
across the junction for simplicity (Fig. |7ja)) In presence of 
a slowly varying potential across the GPNJ, the transverse 
modes with finite wave- vector ky act like waveguides with 
a cut-off frequency and thus need to overcome a real tunnel 
barrier. A smooth GPNJ thus selectively transmits the low 
angle electrons |49|, with suppression of the higher angle 
modes with increasing length of the junction transition re- 
gion. The fundamental normal mode still does not see a bar- 
rier and thus transmits perfectly, reflection once again elim- 
inated by a pseudospin dominated selection rule, retaining 
the Klein tunneling properties of GPNJ. 

The conductance of a symmetric (equal and opposite 
doping on both sides) smooth GPNJ, using the expression 



from Eq. [37] becomes 



4q^ 



I 

j-i 



TikfdO^ 



AG 




(46) 



where Gq = {Tckpd)~^ is the effective critical angle imposed 
by the smooth junction | [49) and AG = Aky / {kcosG) is the 
angular separation between adjacent modes. Numerically cap- 
turing this exponential decrease in conductance with gate 
split distance d requires simulating a long device. The defi- 
nitions of the unit cell Hamiltonian are now changed by the 
following. 

Aka 



A -ika 



(47) 



where is the layer to layer coupling matrix in the trans- 
verse direction. This will describe a layer Hamiltonian with 



periodic boundary condition. Similarly we use a^^ to calcu- 
late gk as described in section II. The Tav vs d is shown in 



inset of Fig. 12 with d extended in the micron regime. 



4 Conclusions - the role of numerical simulation 

Graphene and its heteroj unctions constitute a fascinat- 
ing system that demonstrate rich physics such as anoma- 
lous quantum Hall, chiral tunneling, metamaterial behavior 
and pseudospintronics. Extended to other sub-systems such 
as bilayer graphene, this richness proliferates with Lifschitz 
transitions, anti-Klein tunneling, more anomalous Hall sig- 
natures and so on. The simplicity of the low-energy 2-D 
bandstructures of graphene and its progeny allows pen and 
paper demonstration of these physical concepts. The NEGF 
approach provides a common computational framework to 
verify these concepts atomistically, whereby a small modi- 
fication to the Hamiltonian (e.g. a barrier or a field-induced 
phase in the coupling parameter) can manifest these physical 
concepts with very little effort, even in tougher geometries 
such as multiple junctions and tilted hetero structures. When 
comparing with experiments, however, one needs to worry 
about the sensitivity of these effects to scattering, edge states, 
quantization, smooth junctions, strain, roughness and so on. 
This is where computation plays a definitive role. It was 
hard to anticipate a-priori what the effect of specular edge 
scattering would be on the resistance of a tilted junction, 
but we now know through numerical simulations 1 16] that 
it creates a sweet- spot in the scattering profile that experi- 
ments do not show. Such 'smoking-gun' demonstrations al- 
low us to judge the feasibility of devices that utilize uncon- 
ventional electronic switching in graphene, such as the cre- 
ation of transmission gaps with subthermal switching |T5| . 
We can now atomistically simulate graphene systems whose 
sizes are comparable to experimental dimensions, and di- 
rectly demonstrate many of the physical effects that such 
devices are based on in presence of non-idealities. Further 
extensions will involve the role of self-consistent screening 
(Poisson), charge puddles, defects and incoherent scattering 
from remote optical phonons. 
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